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Abstract. We present an efficient and simple modifica- 
tion of the standard transport algorithm used in explicit 
eulerian fixed polar grid codes, aimed at getting rid of 
the average azimuthal velocity when applying the Courant 
condition. This results in a much larger timestep then the 
usual procedure, and it is particularly well-suited to the 
description of a Keplerian disk where one is traditionally 
limited by the very demanding Courant condition on the 
fast orbital motion at the inner boundary. In this modi- 
fied algorithm, the timestep is limited by the perturbed 
velocity and by the shear arising from the differential ro- 
tation. FARGO stands for "Fast Advection in Rotating 
Gaseous Objects" . The speed-up resulting from the use of 
the FARGO algorithm is problem dependent. In the exam- 
ple presented here, which shows the evolution of a Jupiter 
sized protoplanet embedded in a minimum mass proto- 
planetary nebula, the FARGO algorithm is about an order 
of magnitude faster than a traditional transport scheme, 
with a much smaller numerical diffusivity. 

Key words: Accretion, accretion disks; Hydrodynamics; 
Methods: numerical. 



ratio) decreases as r"'^/^. Since in most cases the "inter- 
esting region" of the grid is located much further than the 
grid inner boundary, the CFL ratio in the region of interest 
is much smaller than unity, which corresponds to a waste 
of computing time, and, as we are going to see below, to an 
enhanced undesirable numerical viscosity. The most obvi- 
ous solution to get rid of such a limitation is to work in the 
comoving frame. Unfortunately, most finite-difference HD 
eulerian codes require an orthogonal system of coordinates 
(Stone & Norman, 1992), which makes them unsuitable if 
one wants to work in the comoving frame in a differen- 
tially rotating disk, and even a non-orthogonal grid eule- 
rian code would be unable to track accurately the fluid 
motion after a few orbits, due to the strong winding of 
the coordinate system. On the other hand, one can adopt 
a Lagrangian description of the disk (Whitehurst, 1995), 
but the implementation is much more tricky and difficult. 
Furthermore, the geometry of an accretion disk provides 
a polar mesh as a natural grid. We describe hereafter a 
simple method which enables one to work on a fixed polar 
grid and to get rid of the CFL condition on the average 
azimuthal velocity at each radius. 



1. Introduction 

We want hereafter to model the hydrodynamical (HD) 
evolution of a disk described on a fixed polar eulerian grid. 
For the sake of simplicity we are only going to deal with 
a two dimensional Keplerian disk, but the algorithm can 
be extended with little additional effort to any gaseous 
thin or thick disk in differential rotation. Usually in this 
kind of numerical simulations the timestep is limited by 
the Courant Friedrich Levy (CFL) condition at the in- 
ner boundary, where the motion is fast and the cells are 
narrow. Indeed, the ratio of the distance swept by the ma- 
terial in one timestep to the cell width must be lower than 
unity over the whole grid, otherwise a numerical instabil- 
ity occurs (i.e. non physical short-wavelength oscillations 
appear, grow exponentially and spoil the model). In a Ke- 
plerian disk this ratio (which we call hereafter the CFL 
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2. Notations and standard method 

We consider a polar grid composed of Ng sectors, each 
one A9 = j^ wide, and Nr rings, with separations at 
radii -Ri(o<i<„^)- The inner boundary is then located at 
the radius Rq, and the outer one at the radius Rn^- 
The density (and the internal energy if needed by the 
equation of state) is centered in the cells, and is de- 
noted {'Sij)(ij)e[o,Nr.-i]x[o,Ns-i]- The radial velocity is de- 



noted 



and is considered centered in azimuth and half- 



centered in radius (applied at radius Ri, i.e. at the inter- 
face between the cells [i, j] and [i — 1, j]). In a similar way, 
the azimuthal velocity is denoted vf,, and is considered 
centered in radius and half-centered in azimuth (i.e. at the 
interface between the cells [i,j] and [i,j — 1] ; throughout 
this paper the algebra on the j coordinate is meant in 
Z/NgZ to account for the periodicity in azimuth). Usu- 
ally in a finite difference code the timestep is split in two 
main parts (Stone & Norman 1992). The first part is com- 
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posed of eulerian substeps which consist in updating the 
HD quantities through the source terms in the evolution 
equations, and which include all the physical processes at 
work: pressure, gravity, viscosity, etc., and which can for- 

E 

mally be described by the transformation f ^ ^", ^ being 
any HD field on the grid. The second part is the transport 
substep, in which the quantities are conservatively moved 
through the grid according to the flow [(«[,■ )°, (wfj)°], and 

which can be formally represented as £,"" -^ £,'' ^ S,'^ , 
where £,~^ denotes any HD field after a whole timestep is 
completed, and R and T denote respectively the radial 
and azimuthal transport operators, which can be alter- 
nated every other timestep. The CFL condition comes 
both from the source part and the transport part, and 
the most stringent restriction is given by the T-substep, 
due to the unperturbed azimuthal flow. Classically, the 
azimuthal transport can be written as: 






Ai_. b,*/vO' 0a 



i. 



b,*/v 
ij + 1 






l) 



(1) 



where Ay^ = .^ '+^ /S.9 is the "mean azimuthal width" 
of a cell. Eq. (nl) expresses the balance of the arbitrary 
conservative quantity ^ in the cell [i, j] by computing the 
difference of its inflow at the [i, j — l]/[i,_;/] interface with 
the velocity wfj* and its outflow at the [i,j + l]/[i,j] in- 
terface with the velocity wf'Yj^. Actually we consider the 

flux of the upwinded interfacial quantity ^^'*/'" , where 
the "*/w^°" operator depends on the numerical method 
(donor cell, van Leer, PPA, see e.g. Stone & Norman 1992) 
and on the velocity field w^°. 

3. New azimuthal transport algorithm 

3.1. Overview 

Let us take as an example the angular momentum conser- 
vation equation: 



dJ_ 
dt~ 



1 d{rv'J) 
r dr 



azim. transport rad. transport 




Source terms (2) 



where J = prv^ . The transport equation of any HD quan- 
tity S, will look the same as the L.H.S. of eq. (||). 

Now without loss of generality we can rewrite eq. (0) 
as: 



a source step, 
a radial transport step, 

an azimuthal transport step with the velocity v^ — w , 
which we are going to call the azimuthal residual ve- 
locity, 

and an additional step which corresponds to the fol- 
lowing PDE: 



dJ_ v^dJ_ 
'dt^Vde ^ 



(4) 



It is an easy matter to check that the solution of this 
last equation can be written in a general way as: 



j{e,t) ^ J 



,0 



(5) 



dJ_ ld\{if_u)J] udJ_ 19(rVV) 
dt r 



89 



89 



8r 



Source terms (3) 



where u can be any quantity which does not depend on 9. 
No assumption has been made on the behavior of J up to 
this point, and eq. (pi) and (ph are strictly equivalent. If we 
take u to be the average azimuthal velocity v^ , then eq. (0) 
can be described as a composition of different steps, and 
each of them can be worked out independently with the 
well-known operator splitting technique: 



which means that the solution of this equation at any 
time t looks like the initial profile (t = 0), except for a 
shift — Jq v^dt/r in azimuth. It should be noted that 
this is true whatever the profile of J, which can even 
contain discontinuities (i.e. shocks). In particular no 
assumption has to he made on the linearity of the flow 
(i.e. on the relative amplitude of the perturbed quan- 
tities). 

A qualitative reason of why such a decomposition is 
valid is that the time evolution of the HD quantities can be 
described either by an observer sitting on a ring of radius 
r which rotates at any instant in time with the average 
azimuthal velocity, or by an observer at rest in an inertial 
frame. Now the time evolution of the system is of course 
observer-independent, which is why their observations are 
reconciled through the simple shift described by eq. (P). 

The idea on which the FARGO algorithm is based on 
is precisely to evolve the HD quantities through opera- 
tors which mimic in a discrete way the different terms of 
eq. (^). The source step, the radial transport step and the 
residual azimuthal velocity transport step are performed 
in a standard way (see e.g. Stone & Norman, 1992). Now 
the last step in the operator-splitting described above, 
which corresponds to a simple shift which amounts to be 
tJ At/r in one timestep, can be implemented in such a 
way that the matter can sweep an arbitrary number of 
cell widths in one timestep. 

In order to lay down the basic mechanism by which 
FARGO works, let us take the following concrete exam- 
ple. We assume that, after the classical substeps (which 
are the source step, the radial transport and the residual 
azimuthal velocity transport), the material at a given ra- 
dius r has to be shifted by 4.7 cells in one timestep (which 
means that lfAt/rA9 = 4.7). What is actually done is 
that 4.7 is decomposed as 4.7 = —0.3 + 5, i.e. the nearest 
integer and a remainder which by construction is lower or 
equal to 0.5 in absolute value. In the first substep of this 
shift step the material is shifted by this remainder (here 
—0.3), which can be achieved through a classical trans- 
port method since the remainder is lower or equal to 0.5 



F. Masset: FARGO : a fast eulerian transport algorithm for differentially rotating disks 



in absolute value (it has to be < 1 in order for the stan- 
dard transport method to be possible), with the additional 
simplicity that the corresponding velocity field is uniform 
(which is actually why shift and transport happen to co- 
incide in this special case, since there is no compression 
in the corresponding flow). The second substep just cor- 
responds to an integer number of cells shift, which is done 
in our example simply by copying the content of cell j into 
cell J -t- 5, for any j. 

A more formal and detailed description of the FARGO 
algorithm is given in the next section. 

3. 2. Mathematical formulation of each step of the 
FARGO algorithm 

In the modified algorithm, the azimuthal transport sub- 
step is split in several parts. We assume that the timestep 
At has already be chosen, and defer discussion of the 
timestep constraints until section [3.3| . We first compute 
the average azimuthal velocity at each radius: 



1 



ea 



Ns-1 
3=0 



(6) 



We then introduce the residual velocity: ff,'°^ — vfj^ — v^ 
and the "shift number" at each radius : 



J = E 


' fl At ■ 




[ ' AyJ 



(7) 



where E[X] denotes the nearest integer to the real X. We 
define the constant residual velocity to be: 



Vi = V, - 71,: 



Ay. 
At 



Hence the total velocity can be expressed as: 



v'^°- = v^^H -I- V^ 



(8) 



(9) 



where the "shift velocity" v\^^ — n.i-^ corresponds to a 
uniform shift of rii cells over one timestep. 

We first transport the HD quantities according to the 






ij + 1) 



(10) 



then to the uniform flow v 






,,ecr 






Ay. ^^^^" 



e 



ij+1 



') 



We split the first part of the transport into two parts {v^'-''^^ 
and v^'^'^) instead of using a single transport step with the 



velocity v' 



Ores 



,,Scr 



in order to ensure (as can be checked 



below given the timestep constraints) that in each of these 
transport substeps the material sweeps at most half a cell 
(it could sweep up to one cell, but for reasons which will 
become clear in section |[ we prefer to take a half cell lim- 
itation), and in order for the continuity considerations of 



section 3.4 to apply. Finally, the quantities are transported 



along the ii^^^ uniform flow 






(12) 



Only the first two parts of this transport step introduce 
some numerical diffusion. The last one, given by eq. (|l2|), 
which in many cases corresponds to the largest part of 
the motion, does not introduce any numerical error, since 
it just corresponds to a circular permutation of the grid 
cells, or in other words it is just an integer discrete version 
of the shift given by eq. (||) . 

A precise quantification of the lower numerical diffu- 
sivity of FARGO is beyond the scope of this paper though. 
An extremely rough estimation can be done in the case of 
the comparison of a standard method (in which the effec- 
tive CFL ratio is a sizable fraction of one) and a FARGO 
method for which n. 7^ 0. If we assume that numerical ef- 
fects will behave in azimuth as a physical viscosity would 
do, then the effective numerical viscosity in FARGO is 
about ni/C'o times lower than the standard method's one, 
where Co is the CFL standard dimensionlcss limitation 
factor, which is detailed in the next section. Neverthe- 
less a variety of numerical experiments can be found be- 
low which all show that FARGO 's numerical diffusivity is 
smaller than the standard method's. 



3.3. Tim,estep limitation 

In the standard transport method, the timestep limita- 
tion arises from the combination of four different con- 
straints (see e.g. Stone & Norman 1992), namely the fact 
that a flow advected test particle in cell [i,j] should not 
sweep a distance longer than Ay, in azimuth nor longer 
than i?i-fi — Ri in radius over one timestep (which in- 
troduces the limit timestep St2 and St^ in Stone & Nor- 
man's paper), and that the wavefront of any wave present 
in the system should not travel across a whole cell over 
one timestep (Richtmyer & Morton, 1957), which corre- 
sponds to the limit timestep Sti in Stone & Norman's pa- 
per. The last constraint comes from a stability limit aris- 
ing from the viscosity (numerical or physical). With the 
modified azimuthal transport algorithm, the constraint on 
the azimuthal motion has to be modified slightly. Fol- 
lowing Stone & Norman's notation, instead of writing 



(11) St'^ = Ayi/vf"-, we write 



5t'i 



Ay, 



Ay. 



,,ea 



(13) 



which means that the timestep limitation comes now from 
the perturbed azimuthal velocity, which results in a much 
higher absolute value of St^. Another limitation arises 
from the shear. Indeed we do not want the shear to dis- 
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connect the two neighboring cells [i,j] and [i + l,j] after 
one timestep. We write this condition as : 






(14) 



Following Stone & Norman's notations, we finally adopt : 



mr' + {5fi^j''Y"] 



3.4- Continuity 



At each timestep, Nr values of rii (with i e [0,Nr — 1]), 
used in eq. (|l^), are computed using eq. @). These inte- 

3/2 

ger values scale roughly as i?j . The shift on the central 
parts generally amounts to several cells over one timestep, 
while in the outer parts rii is small, and possibly zero. 
One can wonder whether or not problems may arise at 
the radii Ri where rii 7^ ?^i-i (i-e. at radii where the az- 
imuthal shift corresponding to the third substep of the 
transport step is discontinuous). More generally we want 
to examine the question of the continuity of ^^^ with re- 
spect to vfAt. In order to check for this continuity, we 
assume v^ = (iV + 5 + e) -^, where N is an integer, and 
we work out the behavior of G (e) in the vicinity of e = 0. 
Since we have to use the explicit form of the "*/t;^°" op- 
erator, we adopt the van Leer algorithm (van Leer, 1977), 
which is widely used. Some straightforward algebra leads 
to : 






ij-N-1 



2 *^ ) i^tj-N 



i 



ij-N-1 



- - -e 



~^ydiij-N 



d^?- 



ij-N-l) 



(16) 



both for e > and e < provided |e| < -^ and where the 
operator "d^" is the van Leer slope. Eq. ( p^ shows that 
the field $^f, is a continuous function of e and hence of ijf . 
In particular no special problem is to be expected from 
the discontinuities of rii across the disk. 



3.5. Operators swapping 

As we said in section g, it is a common practice to al- 
ternate the radial R and azimuthal T transport operators 
every other timestep. In this modified algorithm, R should 
usually be applied first, unless the velocity field is updated 
just after applying the T operator from the new momenta 
and new density fields, or unless special care is devoted 
to the j indices. Indeed swapping blindly the R and T 
operators would result in moving radially the matter with 
the radial velocity it actually has ^ Ui cells upwards, and 
would quickly end in a non-physical staggering everywhere 



4. Mono-dimensional tests 

4-.1. General considerations 

In order to validate this modified transport algorithm, we 
present some ID tests, and we compare the results of the 
standard method and of the FARGO method on a realistic 
test problem. We solve simultaneously the continuity and 
Navier Stokes equation for an isothermal gas (which has 
a non- vanishing but small kinematic viscosity) : 



(15) dp d{pv) 

dt dx 

dv dv 

dt dx 







dp 
dx 



d'^v 
dx"^ 



(17) 



(18) 



We assume that at rest the system has a uniform den- 
sity po and sound speed Cs. The waves which can propa- 
gate in this system have the following dispersion relation- 
ship: 



iO 



±\ k^cl 



kh 



k-v 



or: Lo = ±kcs — i—— n v <^ v\nji = —r- (19) 

2 k 

which reduces to the standard dispersion relation for 
an undamped acoustic wave u = ±kcs provided the sys- 
tem is evolved for a time small compared to the damping 
timescale r = ^ . This will be the case for the results we 
are going to present below, so that any apparent damping 
of the waves has a numerical origin. We do the following: 

1. We first analyze the propagation of a sound wave in 
the matter frame, i.e. we take as initial conditions: 

p{x) = spo cos{kx) and v(x) = scg cos{kx) (20) 

where s is the wave relative amplitude. The polariza- 
tion adopted corresponds to a rightwards propagating 
wave. According to eq. (|l9|) , it propagates with a phase 
velocity which is 5R (^) = c^. We study this propaga- 
tion with the standard transport algorithm (we are in 
the matter frame so there is no systematic average x- 
velocity, hence no need for a FARGO algorithm). We 
check that in this case the solution we get is accurate 
by varying the timestep and checking that the solution 
has converged. 

2. We then turn to a case where the setup is slightly mod- 
ified. We take: 

p{x) = spo cos{kx) and v{x) = vq + scs cos(A:a::)(21) 

where vo is a constant, which we choose much bigger 
than Cs (which would correspond to the conditions of 
a thin keplerian disk, for example). The evolution of 
the system from this setup ought to be the same as 
before, since it merely corresponds to the same phys- 
ical situation, but described from a frame moving at 
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a constant speed —vq w.r.t the first one, so one can 
invoke Galilean invariance to conclude that the wave 
profile evolution has to be the same. So any "good" 
algorithm should approach as closely as possible the 
results of the matter frame simulations. We show that 
this is not quite the case with the standard transport 
method, which suffers from quite a high numerical dis- 
sipation, whereas FARGO behaves much better (not 
to mention its much faster execution). As a side result 
we also show that in this problem taking a CFL effec- 
tive ratio (for the standard transport method) bigger 
than ^ leads to an artificial and non-linear increase of 
the wave profile, and hence has to be avoided. 



4-. 2. ID Numerical results 

We deal with a ID grid composed of Ng = 200 cells, with 
periodic boundary conditions. The cell width is Ax = 
0.0314, the isothermal sound speed is Cg — 0.04. The equi- 
librium density is Sq = 6 • 10~^. These parameters corre- 
spond roughly to the ones used in the numerical study of 
a protoplanet on a circular orbit at 5 A.U. embedded in a 
minimum mass protoplanetary disk (Hayashi et al. 1985 or 
Bryden et al. 1998), that are described in section ^, when 
the central star mass and the protoplanet orbit radius are 
taken to be respectively the units of mass and distance. 
We present the results of different test runs in fig. |l|. The 
thick solid line represents the initial profile, which cor- 
responds to a rightward propagating acoustic wave, with 
wavelength A = 40Aa:: = 1.256. The relative amplitude of 
this sound wave is s = 10^^. The thick dashed line rep- 
resents the density profile at time to — 220, i.e. after the 
wave has traveled Cgto/X = 7 times its own wavelength, 
when studied in the matter frame, i.e. when the velocity 
at t = is set to be only the perturbed velocity associated 
to the sound wave. The thick dashed profile is obtained 
with the standard transport algorithm (there is no need 
for the modified one in this case since we work in the 
matter frame), with a timestep At — 5 • 10~^. The curves 
obtained by choosing a much smaller timestep appear to 
coincide exactly with this one, hence we can consider this 
thick dashed line as the actual state the system must have 
at the date to- This profile does not exactly coincide with 
the initial one because to is ^ i of the profile steepening 
time tjjs ~ T^- 

y^ 2csS 

Now if we just change the initial velocity by uniformly 
adding 1.0 to them at t = 0, which means that we are no 
more in the matter frame, and we still work with the stan- 
dard transport algorithm, then we get the dotted profile, 
which has '^ 1/5 the amplitude obtained from the com- 
putation in the matter frame. In this run the CFL ratio 
is vAt/Ax = 0.16. In order to check the timestep de- 
pendency of this result, we redo this test with twice as 
smaller a timestep {At = 2.5 • lO"'^) and we get the dash- 
dotted profile, which has about twice as smaller a density 
contrast than the previous curve. Note that if this effect 
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Fig. 1. Compared evolution of an acoustic wave evolved 
with the standard transport algorithm and with the mod- 
ified transport algorithm. We plot only two of the five 
wavelengths, i.e. 80 cells out of 200. Due to numerical ef- 
fects the phase velocity of all these profiles do not exactly 
coincide with c^, so that after a time to their phases do not 
coincide . For this reason the profiles have been shifted so 
that they have all approximately the same phase in order 
to improve the clarity of the plot. 



were to be due to a physical kinematic viscosity v, then 



its value should be : v 



X^ log 5 
2T^Ha 



5.8 • 10 , much higher 



than the expected viscosity in a minimum mass protoplan- 
etary disk {v ^ 10~^ in our dimensionless units). Now, 
instead of decreasing the timestep, we increase it and set 
At = 2.0 • 10-2 (hence the CFL ratio is about 0.64). We 
then get at time to the dot-dot-dot-dashed profile, which is 
not numerically damped but slightly amplified. With such 
a large timestep, we can use the modified transport algo- 
rithm, which in that case corresponds to a rightwards one 
cell shift and a leftwards normal transport with a remain- 
ing CFL ratio of 1 — 0.64 = 0.36. In that case we get the 
thin long-dashed profile. If we use the modified FARGO 
transport algorithm, we can still increase the timestep. 
The thin solid profile and the thin short-dashed profile 
have been obtained respectively with At = 4 • 10"^ (ef- 
fective CFL ratio ~ 1.3) and At = 1.2 • 10"^ (effective 
CFL ratio ^ 3.8). We clearly see from these results that 
the FARGO transport algorithm leads to less numerical 
dissipation than the standard transport. From the first 
two tests in the non-comoving frame, one can conclude 
that increasing the number of timesteps over a given time 
interval with the standard transport algorithm increases 
the numerical dissipation (if the grid is moving w.r.t the 
matter frame with a velocity uq 7^ and if the main part 
of the velocity comes from ug). A simple explanation for 
the lower numerical dissipation of the FARGO algorithm 
is that it requires less iterations as the timestep increases, 
and since most of the distance swept is achieved through 
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an exact shift (a circular permutation) , the numerical dis- 
sipation has to decrease as the timestep increases. 



5. Two-dimensional example : the embedded 
protoplanet problem 

We show in this section the validity of the modified trans- 
port algorithm when applied to the interaction of a Jupiter 
sized protoplanet with a minimum mass protoplanetary 
disk in which it is embedded. The perturbed potential as- 
sociated with the planet excites spiral density waves in the 
disk, which propagate away both inwards and outwards, 
with a pattern frequency equal to the planet orbital fre- 
quency. The spiral waves interact with the disk and give it 
the angular momentum they removed from the planet, and 
eventually open a gap centered on the planet orbit, pro- 
vided the planet mass is high enough (Papaloizou & Lin, 
1984). We present a run with a one solar mass primary, 
one Jupiter mass protoplanet initially on a fixed circular 
orbit at ro = 5 A.U. embedded in a standard protoplane- 
tary nebula whose parameters have been mentioned above. 
The grid has an inner radius at 2 A.U. and an outer radius 
at 12.5 A.U. The sequence (-Ri)ig[o^jVr] is equally spaced, 
with Nr = 49 ; The grid has iV^ = 143 sectors, it is fixed 
in a non-Galilean non-rotating frame centered on the pri- 
mary. Its outer boundary is rigid and its inner boundary 
allows outflow but no inflow. The disk aspect ratio is set 
to 4 • 10^^ everywhere. The planet perturbed potential is 
smoothed on a length scale which amounts to 40 % of the 
Roche radius. In eq. ( [14|) we choose Co — 0.5. We plot on 

fig. y the quantity e^ = ^ after 2.86 orbits. This quan- 
tity represents the effective CFL ratio. With the standard 
transport algorithm this ratio is bounded by Cq. 



algorithm in this case results in a speed-up by a factor ~ 8 
of the computation. One can note that the difference in e^ 
between the innermost ring and its immediate neighbor is 





Fig. 2. Number of cells crossed during one timestep. See 
text for parameters. The inner boundary is at the left 
(high values) and the outer boundary at the right (low 
values). 



We see that the innermost ring sweeps almost four cells 
on one timestep, hence the use of the FARGO transport 



0.5, which is the maximum allowed by eq. (14). Indeed the 



timestep in this run is shear-limited, and the constraint on 
the residual velocities only would lead to an even bigger 
timestep, since as one can see the residuals of the distance 
swept over one timestep amounts to far less than 1/2, even 
in the vicinity of the planet. Indeed, runs performed with 
a logarithmic polar grid {i.e. with Ri^i/Ri constant), 
which have a smaller value of Ri+i — Ri in the inner part, 
have shown to allow a speed-up by a factor ~ 30 w.r.t. 
the standard method. 

In order to see how numerical viscosity affects the disk 
response in both cases, we plot on the fig. || the disk den- 
sity after 28.6 orbits, obtained from different algorithms. 
The left plot corresponds to a non-rotating frame stan- 
dard transport run, while the middle plot represents a non- 
rotating frame FARGO transport run, and the right plot 
represents a standard transport run in a frame corotating 
with the planet (hence the planet is fixed with respect to 
the grid, so we expect from the results of section ^ the 
density response in the vicinity of the planet to be given 
with a high accuracy). Note that special care has to be 
devoted to the treatment of the Coriolis force in that case 
in order to conserve exactly the angular momentum and 
then to avoid a spurious outwards transport in the disk 
(Kley, 1998). We clearly see that the global spiral pattern 
excited by the protoplanet in the disk is identical in the 
three cases, though the response in the immediate vicinity 
of the planet is much more spread out in the non-rotating 
frame standard accretion case (left plot), and that the 
most sharply peaked response is achieved through the use 
of a corotating frame (right plot), as expected. Indeed, we 
plot on fig. H a cut of the disk density at the planet radius 
in the three cases. The solid line represents the FARGO 
transport result, and the dot-dashed line the corotating 
frame result. They both have the same width, though the 
maximum of the density in the corotating case is higher. 
The dashed curve represents the result of the standard 
transport in a non-rotating frame. Its width is about twice 
as large as the other curves' width, and we also see that 
numerical effects in that case lead to additional leading 
and trailing material (near cells number 65 and 77), and 
to a smaller density peak value. 

The FARGO plot on figure || exhibits at its inner 
boundary an oscillatory behavior which originates from 
three combined effects. First, this is a shear-limited run 
— see eq. (p^ and fig. ^ — ; if we change the 0.5 fac- 
tor in eq. (|14|) to 0.3, this oscillatory behavior disappears 
(hence in any high resolution run, where the algorithm 
is most likely to be residual velocity limited rather than 
shear-limited, it never turns up). Second, the inner grid 
has strongly radially elongated cells. If we take a log-grid 
(see e.g. Nelson, 1999, or the example below), where the 
cells are almost "square" everywhere, this behavior is not 
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observed, even if the run remains shear-limited. And fi- 
nally, we have a steep density gradient close to the inner 
boundary. If the inner boundary was closed and hence if 
we had no density gradient, this oscillatory behavior would 
never appear. In all the cases where it was observed this 
behavior always disappeared after a few tens of dynamical 
times. 

It should be noted that the numerical damping ob- 
served in the non-stationery frame in section ^ occurs 
both in the non-rotating frame and corotating frame (far 
from the coorbital region) standard method runs. Hence 
the amplitude of the protoplanet triggered density wave is 
marginally higher in a FARGO run at the inner boundary. 
Both this reason and the effect we noticed in the previous 
paragraph lead to a marginally higher mass loss through 
the inner boundary, at least during the first stages of the 
evolution of the system, which results in the darker band 
at the inner boundary in the middle panel of figure ^. 

We present in figure ^ the results of three runs (non- 
rotating standard and FARGO, and corotating standard), 
which describe the same physical system as before af- 
ter the same amount of time, but with a grid for which 
Nr = 70, Ns = 180, Rrmn = 0.25 and Rmax = 2.5, and 
with a geometric sequence for {Ri)i£[o^Nr-] (hence it is a 
log-grid, and everywhere its cells are almost "square"). 
One can check on these plots that there is no oscillatory 
behavior in the FARGO results (this time the cells are no 
radially elongated near the inner boundary), whereas the 
run is still shear-limited. Furthermore, as stated above, a 
careful look at the inner spiral structure shows that it has 
a slightly higher amplitude in the FARGO case. 




Fig. 3. Disk density Ey ; j is in abscissa and i in ordinate. 
The left plot has been obtained by a non-rotating frame 
standard method, the middle one by a non-rotating frame 
FARGO transport method and the right one by a coro- 
tating frame standard method. Since each of these plots 
is approximately square, any circular feature in the disk 
should appear on the plots as a 1 : 3 vertical ellipse. This 
is not quite the case of the material surrounding the planet 
in the left panel, which leads to the conclusion that in a 
non-rotating frame standard transport method, the mat- 
ter is artificially elongated along the orbital motion. The 
FARGO case, in the middle panel, shows much better be- 
havior, and the coorbital material has a distribution which 
looks very much like the right panel one. 




Fig. 4. Disk density Sy ; j is in abscissa and i in ordinate, 
for the log-grid runs described in the text. The left plot has 
been obtained by a non-rotating frame standard method, 
the middle one by a non-rotating frame FARGO transport 
method and the right one by a corotating frame standard 
method. The same comments as in figure |3| apply here. On 
this specific example, the FARGO run turned out to be 
17 times faster than the standard run in the non-rotating 
frame, and 15 times faster than the standard run in the 
co-rotating frame. 
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Fig. 5. Disk density cuts at the planet radius. The solid 
line represents the FARGO transport case, the dashed line 
represents the standard case, and the dot-dashed line rep- 
resents the corotating frame result. The dotted line indi- 
cates the unperturbed surface density. Note that the local 
maxima at j ~ 46 and j ~ 94 correspond to a tempo- 
rary residual accumulation of material at the L4 and L5 
Lagrange points of the protoplanet. 



From the results depicted in figures 0, and 0, one can 
deduce that the FARGO transport algorithm on this par- 
ticular problem is much closer than the usual standard 
transport algorithm to the exact solution (which must 
closely resemble the results given by the corotating frame 
run, at least in the coorbital region, since in section H we 
have seen that one needs to be in the comoving frame in 
order to get accurate results even in the limit of a van- 
ishing timestep) . Another quantitative evaluation of the 
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FARGO algorithm consists in monitoring the accretion 
rate onto the planet as a function of time. We present on 
figure H the accretion rate onto a one Jupiter mass pro- 
toplanet embedded in a minimum mass protostellar disk 
with no initial gap. The disk parameters are the same as 
before, as well as the grid resolution (arithmetic radial 
spacing with Nr — 49 and Ng = 143). Three runs are pre- 
sented with three different schemes: the standard method 
in the rotating frame, which gives, according to section 0, 
the most accurate results, the standard method in the 
non-rotating frame, and the FARGO method in the non- 
rotating frame. We use a slightly different accretion proce- 
dure than the one described by Kley (1999). We see from 
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Fig. 6. Accretion rate as a function of time onto a one 
Jupiter mass protoplanet with three different methods. 
See text for details. 



the curve obtained in the corotating frame that the accre- 
tion rate is about 1.6 • 10"'^ Mj. orbit" after 400 orbits. 
This is in relatively good agreement with Kley's results, 
who gets slightly more than 2.0-10"^ Mj.orbit"^ after 400 
orbits in a similar run, but with a different grid resolution 
and a slightly different accretion protocol. We see from 
figure rathat the accretion rate in the non-rotating frame, 
with a standard method, is smaller than in the rotating 
frame run, by a factor ~ 2. The fact that the accretion 
rate is slower in this case was to be expected from the 
curves of fig. ra. Now the run with the FARGO algorithm 
leads to an accretion rate which is between the rotating 
frame results and the non-rotating frame standard trans- 
port results, and which are closer to the rotating frame 
results. From these considerations again we see that the 
FARGO transport leads to a smaller error w.r.t. the ro- 
tating frame results. The point here is that the FARGO 
transport algorithm is about one order of magnitude or 
more faster than the corotating frame standard transport 
run, and that the corotating frame is suitable only to the 
study of a protoplanet on a fixed circular orbit. From these 



remarks it clearly appears that the FARGO transport al- 
gorithm is particularly well suited to the study of the pro- 
toplanet orbit long-term evolution. FARGO has already 
been used to study the migration and mass accretion of 
a Jupiter sized protoplanet in a protoplanetary disk. It 
has been extensively tested against existing independent 
codes, which use the standard transport algorithm. It has 
proven to give very similar results, and the slight differ- 
ences which remain between these codes and FARGO can 
all be understood in terms of FARGO 's lower numerical 
diffusivity (Nelson et al. 1999). 

6. Conclusion 



The FARGO algorithm for the azimuthal transport turns 
out to be able to speed up by about an order of magnitude 
the numerical simulation of a differentially rotating disk, 
with a smaller numerical viscosity than the usual trans- 
port algorithm. It has been validated by many tests on the 
embedded protoplanet problem. It is worth mentioning 
that the FARGO transport algorithm must be used with 
a good understanding of the physical processes at work in 
the system. In particular, the timestep given by eq. (15) 



must be short compared to all the physical time scales 
relevant for the system. In the case we have presented 
in this paper this is automatically ensured by the set of 
eq. ( p^ ) to (|l^), but if additional physics is to be added 
(magnetic field, radiative transfer, etc.), the timestep limit 
needs to be carefully worked out. Furthermore, no advan- 
tage is gained in using FARGO in problems where the 
perturbed velocity is comparable to the rest velocity. It 
is the case for instance of the gas flow in a galactic bar. 
This does not mean that the FARGO algorithm leads to 
wrong results in that case, but simply that it will not be 
better than a standard method, both in terms of numer- 
ical diffusivity and execution time. On the other hand, 
the FARGO algorithm appears to be very well suited to 
all the cases where the perturbed velocities in any differ- 
entially rotating disk are small compared to the unper- 
turbed velocities, which does not mean that the problem 
under consideration has to be linear; indeed the relative 
perturbed amplitude can be arbitrarily high (see e.g. sec- 
tion ra in which the protoplanet wake generates shocks in 
the disk). More generally the FARGO algorithm can be 
used to describe the HD evolution of any sheared fluid on 
a fixed orthogonal eulerian grid. 
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